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Grain boundaries (GBs) are often the preferred sites for void nucleation in ductile metals. However, it has 
been observed that all boundaries do not contribute equally to this process. We present a mechanistic 
rationale for the role of GBs in damage nucleation in copper, along with a quantitative map for predicting 
preferred void nucleation at GBs based on molecular dynamics simulations in copper. Simulations show a 
direct correlation between the void nucleation stress and the ability of a grain boundary to plastically deform 
by emitting dislocations, during shock compression. Plastic response of a GB, affects the development of 
stress concentrations believed to be responsible for void nucleation by acting as a dissipation mechanism for 
the applied stress. 

The microstructure of metals and alloys consists of a network of single crystals (called grains) and the 
boundaries between these crystals (called grain boundaries). While microstructure represents the overall 
uniform description of a material, it is the imperfections or defects like grain boundaries, the weak links, 
which dictate some material properties, especially its propensity to fail. In the case of ductile metals, like Cu, Au, or 
Al, it is known that all grain boundaries are not equal in their propensity for either strengthening or weakening a 
material. However, the reason for this selectiveness is not well understood. To develop next generation predictive 
capabilities for dynamic failure in engineering materials and to design future damage-tolerant materials, it is 
imperative to not just qualitatively understand the reasons behind why all grain boundaries are not equal but also 
to quantify this selectiveness in a manner that can be useful for both design and prediction. In this work, we focus 
on developing a qualitative and quantitative understanding of failure at grain boundaries in ductile materials. 

Fracture in ductile materials is typically characterized by void nucleation, growth and coalescence. A broad- 
base fundamental problem in materials science is where, when, and why voids nucleate and grow in materials as 
an initiation of damage and failure. Previous experimental research has shown that microstructural features such 
as grain boundaries, inclusions, vacancies and heterogeneities can act as initial void nucleation sites 1 " 6 . However, 
we cannot accurately predict the preferred locations and stresses for void nucleation. In fact, larger scale simula- 
tions that study void growth and failure tend to begin with material microstructures where voids have been 
nucleated in a stochastic manner, which is not a true representation of void nucleation in real experiments 7 . This 
random nucleation of voids is tied to a lack of systematic and detailed understanding of the critical mechanisms 
and stresses involved in void nucleation. Direct observations of void nucleation and early stage growth are critical 
to understanding these underlying mechanisms. However, these in-situ observations have long remained a grand 
experimental challenge, largely due to the extremely fine spatial (nm) and narrow temporal (fs-ns) scales 
involved. Hence in this work, we utilize molecular- dynamics simulations, which can provide in-situ, time 
resolved insights, to quantify for the//rsf time the relationship between void nucleation stress and grain boundary 
structure. 

A specific type of fracture in solids, induced by shock compression-and-release termed "spall" motivates this 
investigation 2,8 . During spall, compression waves are generated in both the target and flyer. These compression 
waves are reflected from the target and flyer free surfaces and intersect each other to form a region of tension, 
which can lead to failure 9 . Although, there are differences in stress states, rates of work hardening and operative 
deformation mechanisms preceding damage evolution between various loading conditions, we argue that the 
observations from the current work on dynamic loading provide previously missing insights into the rationale for 
preferred void nucleation sites. 

Specifically, this MD study is motivated by a plethora of experimental work with similar observations regarding 
void nucleation under both dynamic and quasti-static loading 10 " 16 . Dynamic damage in high-purity Cu demon- 
strates, that voids preferentially nucleate at grain boundaries 10 , and specifically certain special boundary types are 
consistently more resistant to void nucleation 11 . Furthermore, it has been shown that boundary structure and not 
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Figure 1 | GB structures at 0 K calculated using an EAM potential for Cu 33 for the boundaries utilized in this work. These structures are viewed along tilt 
axis. The atoms are colored by the centro symmetry 29 parameter where blue represents atoms close to an FCC environment and green represents atoms no 
longer in FCC configurations. The range of the centrosymmetry parameter in the E3 and E5 boundaries is 0-3 and 0-13 in the £11 boundaries. 



just the crystallographic orientation on either side of the boundary is 
important with regard to the early stages of damage 17 . This obser- 
vation is also supported by the work of Wayne et al. 12 who found that 
grain boundaries with certain misorientations, in polycrystalline Cu, 
serve as preferred locations for intergranular damage. These insights 
regarding damage nucleation at high strain-rates apply equally to 
damage at lower strain-rates. Work by Mikhailovskij 13 on bcc tung- 
sten shows that under uniaxial stress conditions, special boundaries 
such as El, E3, and E9 possess higher resistance to failure in com- 
parison with other types of CSL and non-CSL boundaries. 
Experiments by Lim 14 on low-cycle fatigue of polycrystalline fee 
nickel samples showed that low- order CSL boundaries such as E3 
and E5 did not fail during the deformation process. Evrard 15 and 
McMurtrey 16 made similar observations through both experiment 
and simulations showing that special boundaries were more resistant 
to crack nucleation in pre-irradiated austentic stainless steels. These 
results collectively suggest that all grain boundaries, regardless of 
materials or loading condition, are not equal in terms of their pro- 
pensity for serving as void nucleation sites. 

Similar numerical results have been ascertained by molecu- 
lar-dynamics (MD) simulations. Recent MD simulations by Luo 
et al. 18,19 on a £3 asymmetric tilt grain boundary support the experi- 
mental findings listed above by showing that crystal anisotropy and 
grain boundaries both change the response of a material under shock 
compression, in comparison to the response a single crystal, by redu- 
cing both its spall strength and flow stress. Work by Han et al 20 on E3 
coherent and incoherent twin GBs in copper shows that, under 
dynamic loading conditions, void nucleation occurs at the coherent 
twin boundary and not at the incoherent boundary emphasizing yet 
again that not all GBs are equal in terms of void nucleation. 

In previous MD simulation work, we have attempted to under- 
stand and predict void nucleation at grain boundaries by using aver- 
age GB properties 21 " 23 . To help frame the relationship between these 
properties and void nucleation 21 , we considered a standard model for 
the fracture toughness of a material 24,25 : 

Yf = 2Ys-YGB + Y P > (1) 

where y^ y s , y GB and y p are the fracture energy, the surface energy, the 
grain boundary energy, and the energy associated with plastic work, 



respectively. We have examined the importance of y s , y GB , and 
related average properties, such as excess volume to predict the stress 
required to void nucleation at a grain boundary in a ductile metal. 
Those results suggest that there is no direct correlation between grain 
boundary energy 21 , excess volume and the void nucleation stress for a 
ductile material. In a brittle material, where y p = 0, the susceptibility 
of a boundary to fracture can be predicted using the grain boundary 
and the free surface energies. However, the plastic energy term, y p can 
be dominant in ductile materials 10 " 12 . In this work, we determine if 
the plastic work term, y p is a better determinant of the void nuc- 
leation stress for an interface in ductile bi-crystals. 

This paper presents a correlation of the susceptibility of a bound- 
ary to nucleate a void under tensile unloading to its ability to undergo 
plastic deformation under shock compression, using MD simulations 
of copper bi-crystals. Specifically, we simulate the impact of a dissip- 
ative process such as plastic deformation by dislocation emission on 
the development of the stress concentrations believed to drive void 
nucleation. The effort outlined above is the first step towards exam- 
ining and predicting void nucleation at grain boundaries under high 
strain-rate, uniaxial strain loading. This paper presents the simu- 
lation methodologies, the substance of the simulation results, and 
finally conclusions concerning the effect of dissipative plastic 
deformation on the stress to initiate void nucleation at a boundary. 

Results 

GB structures are OK. Four boundaries in Cu were chosen for this 
study - representing a range of GBs encountered in "real" 
polycrystalline materials: 

1. Ell <101> {545J-181} 50.48° asymmetric tilt 

2. Ell <101> {545}-{181} 50.48° asymmetric tilt boundary 
whose structure was disordered by annealing and subsequently 
quenching the Ell asymmetric tilt boundary 

3. E5 <100> {310} 36.9° symmetric tilt 

4. E3 <110> {112} 36.9° symmetric tilt 

The relaxed structures for these boundaries at 0 K are shown in 
Fig. 1 and compare well with other GB structures found in other 
studies 26 " 28 . The atoms in Fig. 1 are colored by the centrosymmetry 
parameter 29 where blue sphere denotes an atom in an environment 
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corresponding to the fee structure and the red sphere indicates an 
environment that is furthest away the from fee structure. The 23 and 
both the El 1 boundaries have pre-existing stacking faults protruding 
from the grain boundary into one of the grains. These stacking faults 
are connected with Shockley partial dislocations that are generated 
due to the dissociation of grain boundary dislocations that reduce the 
stress concentrations at the boundary grain boundary. The average 
calculated properties for these boundaries are available in Ref. 30. 

GB Plasticity under Mechanical Loading. Under shock compres- 
sion, dislocation emission is observed from grain boundaries in the 
perfect bi-crystals. The simulations cells are designed such that other 
complexities like the presence of pre-existing dislocations that can 
lead to pile-ups, triple points and point defects are avoided. Hence, in 
this work, as a starting point, we investigate one particular dissipative 
process, plastic deformation by dislocation emission at a grain 
boundary. 

The extent of plastic deformation by the emission of Shockley 
partials varies in each boundary. The £3 boundary emits regularly 
spaced Shockley partials in both the boundaries as shown in Fig. 2. 
This difference in plastic response is attributed to grain boundary 



structures. This is best highlighted by the 2 1 1 asymmetric tilt bound- 
aries, in Fig. 2, that have different local boundary structures, but the 
same average grain orientations across the boundary, and very dif- 
ferent emission responses. 

We quantify the difference in the plastic response of these four 
boundaries by using the total excess energy (as defined in the 
Methods Section). Figure 3 shows the calculated total excess energy 
(Eq. 2) as a function of time. Before the shock compression wave 
reaches the GB, this excess energy is equal to the grain boundary free 
energy. The arrival times of the shock compression wave at the GB 
vary by 1-2 picoseconds at each boundary because of slight differ- 
ences in the sound speed, in each bi-crystal, due to the varying elastic 
constants (as shown in the inset of Fig. 3). 

As the shock wave reaches the boundary (around 7 ps), the excess 
energy increases sharply. The magnitude of this increase is depend- 
ent on the amount of plastic deformation at a grain boundary under 
shock compression. The plastic work continues to increase steadily 
after the compression wave is reflected from the free surfaces. Even 
though the reflected wave annihilates some of the dislocations emit- 
ted under the compression wave, it also generates more plasticity by 
emitting a different set of dislocations from the boundary. 
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211 asymmetric 
25 symmetric 
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Figure 3 | Total excess energy increases as a function of time after interacting with the shock wave due to plastic deformation. The red line with squares 
corresponds to the 23 symmetric tilt, green with triangles to the 211 asymmetric ordered tilt, blue with circles to 25 symmetric tilt and black with 
diamonds to 21 1 asymmetric disordered tilt boundaries, respectively. The inset shows the slightly different arrival times for the shock compression wave 
at each boundary. 
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Hence, the peak in the excess energy does not necessarily corre- 
spond to the reflection of the compression wave from the free sur- 
faces. Eventually, the annihilation of plasticity is higher than the 
generation of plasticity and results in a drop in the excess energy. 
The increase in the excess energy at later times (after 25 ps) is assoc- 
iated with the nucleation of voids at the grain boundary. However, 
the excess energy at the boundary, immediately prior to void nuc- 
leation is more relevant to understanding and predicting its suscept- 
ibility to void nucleation. 

We extract the times associated with void nucleation to be 
approximately 19.6, 20.2, 21.1 and 19.9 ps for the Ell ordered, 
£11 disordered, Z3 and S5 boundaries, respectively, from our MD 
simulations in which the cells shown in Fig. 1 were loaded to failure. 
The excess energy is then calculated 2-3 ps before voids nucleate at 
the boundaries and is plotted along with the void nucleation stress in 
Fig. 4. This figure shows that there is a direct correlation between the 
total excess energy and the void nucleation stress. It is important to 
note that this is the first time a quantitative relationship has been 
developed and demonstrated between plasticity and the void nuc- 
leation stress. 

Figure 4 shows that grain boundaries that undergo higher plastic 
deformation under shock compression require higher stresses for 
void nucleation. For example, there is a 1.34 GPa difference in the 
void nucleation stress for the two-Ell asymmetric tilt boundaries. 
This difference can be attributed to the fact that the disordered Ell 
boundary can dissipate higher amounts of energy associated with the 
applied stress, by emitting more dislocations as compared to the 
ordered £11 boundary (refer Fig. 4) even though the only difference 
between the two is boundary structure; orientation relationship 
across these two is the same. These results suggest that plastic 
deformation by emission of dislocations from a grain boundary, 
controlled, at least in part, by boundary structure, acts as a principal 
dissipation mechanism for the applied stress and can dictate if the 
stress concentrations required for void nucleation at a boundary are 
achieved. 



can occur, these experimental results amongst many others, indir- 
ectly verify the findings from the molecular dynamics simulations. 

While these results suggest plastic deformation as a principal dis- 
sipation mechanism for the four boundaries studied here, this cor- 
relation has not been rigorously tested for other boundaries and 
materials. To fully quantify the operative deformation and eventually 
damage evolution and failure mechanisms, it will be necessary to 
account for potential other competing processes such as existing 
dislocations, in conjunction with the plastic dissipation phenomena 
explored here. Nonetheless the results of this study for the first time 
quantitatively and qualitatively highlight the importance of plasticity 
via dislocation emission as a dissipation mechanism during void 
nucleation. These findings can be utilized to design a set of criteria 
that can be used to predict void nucleation in a deterministic manner 
and improve our damage prediction capabilities. Finally, this is an 
important step towards modelling the complex damage processes 
and the competition between diverse stress accumulation and dis- 
sipation processes in polycrystalline metals. 

Methods 

All of the MD simulations are based on an embedded- atom method (EAM) intera- 
tomic potential model for Cu developed by Mishin et a/ 33 . This model for Cu has 
proven capable of properly representing a wide spectrum of material properties under 
these moderate conditions, specifically in reproducing the experimental shock 
Hugoniot 34 . MD simulations applied a combination of SOLVER 35 for initially relax- 
ing the grain boundaries at 0 K, and LAMMPS 36 for shock-loading simulations. 

GB structures at OK. In order to obtain a relaxed GB structure at 0 K, we create an 
initial simulation block with two grains separated by a grain boundary. The block is 
periodic in the plane of the boundary (x,z directions) and non-periodic perpendicular 
(y-direction) to it. The grains are sandwiched between two slabs parallel to the grain- 
boundary plane. The atoms in each of these slabs are fixed relative to each other and 
can only undergo rigid body motion. The two grains are free to move and undergo 
displacement in all the three directions during relaxation using a viscous drag 
algorithm. This relaxation technique is equivalent to minimization with respect to the 
specific volume. A standard method of molecular dynamics that generated atomic 
trajectories as a function of time by integration of Newtonian equations of motion is 
used to relax its structures 35 ' 37 ' 38 . 



Discussion 

In this work, we demonstrate by atomistic simulations that there is a 
direct correlation between the stresses required for void nucleation 
and the extent of plasticity at GBs. For the particular boundaries 
studied here, we found that the Z3 symmetric tilt boundary was 
the most resistant to void nucleation due to its ability to efficiently 
dissipate stresses associated with mechanical loading via maximum 
dislocation emission. 

These atomistic results are consistent with experimental findings 
on dynamically loaded Cu that show that £3 twin boundaries are 
more resistant to void nucleation 31 . Indirect experimental findings 
show that grain boundaries that are surrounded by "softer" grains, as 
defined by the Taylor factor, tend to not nucleate voids as easily as 
boundaries that are surrounded by "harder" grains 32 . Softer orienta- 
tion grains promote plasticity whereas harder grains retard it. Since 
Taylor factor is a measure of the ease with which plastic deformation 



Mechanical Loading of GBs. For the shock simulations, the bi-crystal cell is divided 
into two regions normal to the boundary plane: flyer and target as shown in Fig. 5. 
The target region is twice as long as the flyer plate such that the spall plane is located at 
the grain boundary. The shock particle velocity is denoted by u p and in this paper 
it was set to 0.5 km s -1 . Each atom in the flyer plate is initialized with a velocity of 2u p , 
in addition to its thermal velocity, in the velocity component along the shock 
direction, causing it to impact the target. The shock direction is perpendicular to the 
boundary plane. Shock simulations use a NVE ensemble and are performed using 
LAMMPS. Periodic boundary conditions are applied along the x and z-axis and the 
nonimpact sides of the flyer plate and target are free surfaces. 

For the bi- crystals the GB is placed approximately in the middle of the target region 
and the shock waves travel from grain I to grain II, respectively, before being reflected 
from the free surfaces. In our shock loading MD simulations, we use a time step of 1 fs 
and total simulation times of 100 ps. Simulations are initialized at 300 K. 

Stress Required for Void Nucleation. To track the shock front and other physical 
properties during the simulations, the simulation cells were divided into bins of length 
0.55 nm - equal to the cut-off distance for the EAM potential - along the shock 
direction. Properties such as particle velocity and stresses were then averaged within 
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Figure 4 | The stress required for void nucleation increases linearly with 
the calculated excess energy. The points in this plot correspond to the four 
boundaries examined in this study. 
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Figure 5 | The schematic shows the simulation cell used during dynamic 
loading. The red lines represent the grain boundary, the blue lines show 
examples of slip planes along which dislocations can be emitted from the 
boundary. y p and y s are plastic dissipation and surface free energy. The free 
surfaces are perpendicular to the GB (not marked). 
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Figure 6 | Averaged stress as a function of bin position at time, t. The red 

rectangle represents the areas of tension. 

each bin to give the longitudinal shock profile. The atoms are coloured by the 
centrosymmetry parameter to characterize local deformation and structure in the 
simulation cell since this parameter can identify dislocation cores, stacking faults and 
Shockley partials 29 . In this work, void nucleation is defined as formation of a void that 
is coupled (no size or location constraints were applied) with the stress relaxation. 
Nucleation is determined by examining both the quantitative stress data (Fig. 6) and 
qualitative images (Fig. 7). 

The void nucleation stress for the boundaries is calculated as the peak stress, from 
the averages stress profiles, immediately before void nucleation, in the shock dir- 
ection. This is a standard method for calculating stresses associated with void nuc- 
leation in shock loading studies 39 ' 40 . Examples of an average stress profile and 
qualitative images showing the GB before and after void nucleation are shown in 
Figs. 6 and 7. 

Excess Energy. In order to quantify the amount of plastic work, excess energy was 
used as a first order approximation. Excess energy was calculated in each bin using the 
following equation: 

■ntotal _ 1 = 1 



where n is the number of bins, E* ur is the total energy of atoms in bin i, Ep erf is the 
energy of atoms in bin i if they were in a perfect FCC configuration, and A GB is area of 
the GB plane. Excess energy would be approximately zero in bins far away from the 
boundary and have a finite value in the bins that include the GB. The total excess 
energy would be equal to GB energy, prior to deformation. To exclude the affects from 
the thermal and mechanical noise, we only use a limited area around the grain 
boundary to calculate excess energy. This limited region is selected based on the 
centrosymmetry parameter and represents the region with the largest increase in the 
centrosymmetry parameter as compared to a reference. The reference 
centrosymmetry parameter is chosen as the average parameter value in bins far away 
from the boundary at time 0 ps. 

It is important to note that the excess energy contains both the elastic and plastic 
energies. However, a simplifying approximation assuming the elastic component of 
the excess energy to be similar for all the boundaries is made in this study. This is a 
valid approximation since the energy per atom in the compressed region of the 
simulation cell, where the material is only undergoing elastic deformation, varies by 
—0.01 eV amongst the boundaries. 




Figure 7 | Snapshots showing the L3 GB before and after void nucleation. 

The atoms are coloured by the centrosymmetry 29 parameter where blue 
represents atoms close to an FCC environment and green represents atoms 
no longer in FCC configurations. 
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